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ABSTRACT 

Using 1000 ray-tracing simulations for a A-dominated cold dark model in ISato et al.l (|2009l ). we 
study the covariance matrix of cosmic shear correlation functions, which is the standard statistics 
used in the previous measurements. The shear correlation function of a particular separation angle 
is affected by Fourier modes over a wide range of multipoles, even beyond a survey area, which 
complicates the analysis of the covariance matrix. To overcome such obstacles we first construct 
Gaussian shear simulations from the 1000 realizations, and then use the Gaussian simulations to 
disentangle the Gaussian covariance contribution to the covariance matrix we measured from the 
original simulations. We found that an analytical formula of Gaussian covariance overestimates the 
covariance amplitudes due to an effect of finite survey area. Furthermore, the clean separation of the 
Gaussian covariance allows to examine the non-Gaussian covariance contributions as a function of 
separation angles and source redshifts. For upcoming surveys with typical source redshifts of z s = 0.6 
and 1.0, the non-Gaussian contribution to the diagonal covariance components at 1 arcminute scales is 
greater than the Gaussian contribution by a factor of 20 and 10, respectively. Predictions based on the 
halo model qualitatively well reproduce the simulation results, however show a sizable disagreement 
in the covariance amplitudes. By combining these simulation results we develop a fitting formula to 
the covariance matrix for a survey with arbitrary area coverage, taking into account effects of the 
finiteness of survey area on the Gaussian covariance. 

Subject headings: cosmology: theory - gravitational lensing - large-scale structure - methods: numerical 



1. INTRODUCTION 

Weak gravitational lensing by intervening large 
scale structure, the so-called cosmic shear, provides a 
powerful probe of dark matter and dark e nergy. Since 
its fi r st detections by v a rious groups dBacon et all 
20001: iKaiser et all [20001 : iVan Waerbeke et al.l l2000l 
Wittman et al. I2000D . substantial progress has 



been made on both theoretic al and observ ational 
sides (e.g. lHamana et all l2003|: iMassev et al.l 
Jarvis et al.ll2006t ISemboloni et all 120061 : iFu et all 



Schrabback etHT feOlO). 
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Weak lensing has the highest potential to con- 
strain properties of dark energy among other cosmo- 
logical observations, such as type la su p ernovae (e.g 
Riess et al.l 119981: iPerlmutter et all 119991: iHicken et al 



2009), baryon acoustic o scillat ions (e.g. lEisenstein et~ 
20051 : lOkumura et al l I2008T). galaxy clusters (e.g. 



Vikhlinin et al.l 120091: iMantz et alll2010D. if the system- 
atic errors are well unde r control ([Albrecht et al.ll2006l 
l2009trJo~udaki et aT1l2009D . The growth rate of mass clus- 
tering can be measured by "lensing tomogra phy" (e.g. 
iHv] 11993 IHutererl 12001 : iTakada fc .Tainl l20f)l which in 
turn provides tight constraints on the equation of state 
of dark energy. For this purpose, a number of ambi- 
tious wide-field surveys have been proposed, such as Sub - 
aru Hyper Suprime-Cam Survey (|Mivazaki et all [2006). 
the Panoramic Survey Telescope & Rapid Response Sys- 
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tern (Pan-STARRS3), the Dark Energy Survey (DES3), 
the Large Synoptic Survey Telescope (LSST0), the Wide- 
Fiel d Infrared Survey Te lescope (WFIRS10), and Eu- 
clid (|Refregier et a!j|2(JT0h . 

To attain the full potential of upcoming lensing sur- 
veys, it is essential to analyze data with appropriate 
statistical methods as well as to use sufficiently accu- 
rate theoretical models for the power spectrum and/or 
two-point co rrelation function and for the covariance 
matrix (e.g. iHikage et al.l 120111 for a recent develop- 
ment of lensing power spectrum measurement method) . 
Most of the useful cosmological information contained in 
the cosmic shear signal lies in small scales that are af- 
fected by nonlinear clustering. Therefore, non-Gaussian 
errors can be significant in weak lensing measure- 
ments as indicated by several studies ([White fc Hull20"00l: 
Coorav & H"ul 120011; ISemboloni et all 120071: iDore et al l 
2009 1 : lEifier et all 120091: ISato et al l 120091: ITakada fc Jainl 
20091 iLu et al.1120101 : iPielorz et allboioflSeo et al.ll2011l ) 
and also may cause biases in the bes t-fit parameters 
(jHartlap et all 120091: llchiki et all 12009ft . Furthermore, 
future high-precision measurements may require ade- 
quate statistical methods, i.e. accurate likelihood func- 
tion, of weak lensing in esti mating cosmological parame- 
ters (jSato et alll201(ll2l)Tl . 



1 http://pan-starrs.ifa.hawaii.edu/public/ 

2 http: / /www. darkenergysurvey.org/ 

3 http://www.lsst.org/ 

4 http://wfirst.gsfc.nasa.gov/ 
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In the first paper of a series of our works (|Sato et al.l 
2009), we studied the non-Gaussian effects on the covari- 
ance matrix of the cosmic shear power spectrum using 
1000 realizations of ray-tracing simulations. In this pa- 
per, we study the non-Gaussian effects on the covariancc 
matrix of the cosmic shear correlation function which is 
the most conventionally used statistical method in the 
previous measurements. 

This paper is organized as follows. In § [2] we briefly 
review the basics of the cosmic shear correlation func- 
tion and its covariance. In §[3] we show the main results, 
and develop the fitting formula to compute the covari- 
ance matrix of cosmic shear correlation functions for a 
survey of arbitrary area. § [4] is devoted to conclusion. 
Throughout the present paper, we adopt the concordance 
ACDM model with matter density f2 m = 0.238, baryon 
density Sib = 0.042, dark energy density SIa = 0.762 
with equation of state parameter w = — 1, spectral in- 
dex n s = 0.958, the variance of the density fluctuation 
in a sphere of radius 8 /i _1 Mpc erg = 0.76, and Hubble 
parameter h = 0.732. These p arameters are consis tent 
with the WMAP 3- year results (jSpergel et alJl2007ft . In 
our ray-tracing simulation, each realization has an area 
of 25 deg 2 . The detailed description of our ray-tracing 
simulations is described in lSato et al.l (|2009l) . 

2. PRELIMINARIES 

2.1. Real-Space Correlation Function and Its 
Covariance Matrix 

In this section we briefly review definitions of the cos- 
mic shear correlation function and its covariance matrix. 

Since the shear field is a spin-2 field, the field at one 
particular point on the sky carries two degrees of free- 
dom. We can thus define different correlation functions 
from the measured shear field. The most conventionally 
used functions are given in terms of the lensing power 
spect rum as (e.g. Sc hneider et al.1 120021 : iMunshi et al.l 
[20081) : 

e+(*)= r i -^p K (i)j (w), (i) 

Jo 27r 

u&) = l°° l -^p«(i)Mm, (2) 

where P K (l) is the convergence power spectrum (see, e.g. 
Bartelmann fc Schneider 2001; Van Waerbck e fc Mellierl 
20031 for the definition), and Jo(x) and Ji{x) are the 
zeroth and fourth order Bessel functions, respectively. 
Observationally £+ (9) can be measured by averaging the 
product of ellipticity components over all the galaxy pairs 
that are separated by the angle 9 (also see Appendix IX)) . 

In reality the measured £+ and £_ are contaminated 
by systematic errors. Hence it is in practice useful to de- 
compose the measured £+ and £ _ to the lensing-induced 
-E-mode (gradient-mode) correlation function or equiva- 
lently the correlation function of the convergence field, 
£, K (9). Another independent B-mode correlation func- 
tion can be used to monitor residual systematic errors. 
Although £ K and contain theoretically equivalent in- 
formation in the absence of systematic errors, which both 
have the same expression in terms of P K (Eq. [IJ , the two 
show a slight difference when measured from a finite-area 
survey or simulations as we will show below in detail. 



The difference is ascribed to the fact that the conver- 
gence is the projected mass density field, while the shear 
field is a quantity arising from the non-local tidal field 
that is affected by the mass distribution outside survey 
area or simulation area. Therefore in this paper we focus 
on £ + , rather than £ K , to study the covariance matrix as 
the shear field is a more direct observable from actual 
data. 

In Fig. [1] we compare the correlation function £ + (9) , 
measured from the shear fi eld in 1000 ray-tracing simula- 
tions (see lSato~ et al. 20091 f° r t ne details of simulations), 
with the analytical predictions. To obtain the analytical 
predictions we first need to compute the lensing power 
spectrum P K (l), which is given as a projection of the mass 
power spectrum weighted with the radial lensing kernel 
along the line of sight based on the Limber's approxima- 
tion: 

P K (l) = [ S dz W GL (z,z s )P s (k = -^;z] , (3) 
Jo V X(z) J 

w here Wql(z , z s ) is the lensing weight function (see Eq. 9 
in iSato et all [2P09) and Ps (k) is the mass power spec- 
trum . We use the HaloFit fitting formula (|Smith et al.l 
2003) to compute the nonlinear mass power spectrum 
P$(k) for the cosmological model we have assumed. In 
Fig. [T] we show the two analytical predictions. One is 
denoted by the solid curve showing the prediction com- 
puted based on the conventional method; the two-point 
correlation function £+(#) is computed by inserting the 
lensing power spectrum, computed from Eq. ([3]), into 
Eq. |T]). The other is denoted by the dashed curve, where 
we properly take into account the fact that N-body sim- 
ulations used in ray-tracing simulations do not contain 
density perturba tions with length scales beyond the sim- 
ulation box (see lSato et al.ll2009l for details of the sim- 
ulations). For this purpose we imposed P$(k;z) — at 
k < fcf and at each lensing redshift z in computing the 
lensing power spectrum (Eq. |3)), where k{ is the funda- 
mental mode of N-body simulation and given in terms of 
simulation box size L as kf = 2it/L. Note that, since our 
ray-tracing simulations are done in a light-cone configu- 
ration along the line of sight, some large-length modes 
are indeed beyond the area of ray-tracing simulations. 

Fig. [D shows that the two analytical predictions are in 
good agreement with the simulation results at small sep- 
arations. However, the prediction (solid curve) including 
all the modes beyond simulation box overestimates the 
simulation results at large separations, 9 Si 10'. On the 
other hand, the dashed curve, which includes only a finite 
box-size effect of the modes, better reproduces the sim- 
ulation results at the large separations, up to 9 ~ 100'. 
Thus these results imply that the two-point correlation 
function of a given separation angle is affected by a wide 
range of Fourier modes due to the non-local integration 
relation between the real- and Fourier-space modes. In 
particular, even if focusing on the correlation functions 
at large separations in the linear regime, we need to prop- 
erly take into account an effect of the density perturba- 
tions beyond a survey area, which cannot be observed. 
It is also worth noting that, if we impose the cutoff on 
multipole rather than the 3D wavenumber, e.g. including 
only the multipole modes at I > ^ cut = 72 correspond- 
ing to the largest angular mode of our simulated area 
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Fig. 1. — The cosmic shear correlation function, £+(0) (see Eq.[T} 
for source redshift z a = 1.0, measured from our 1000 ray-tracing 
simulations. The error bars at each bin are estimated from the 
scatters (lc scatters) in the measured correlation functions of 1000 
simulations. Note that each realization has an area of 25 square 
degrees, therefore the error bars show the sampling variance for the 
area coverage. For comparison the solid and dashed curves show 
the analytic predictions computed using HaloFit. The difference 
between the two curves is that the dashed curve includes only the 
contributions of mass density fluctuations with length scales within 
the box size of simulations (see text for details) . The dashed curve 
is in better agreement with the simulation result over a wide range 
of separation angles. 

25 sq. degrees, the theoretical prediction underestimates 
too much the simulation result at large separation angles. 

The results in Fig. Q] ca n be compared w ith the re- 
sult of our previous paper (jSato et al.l I2009T ). where we 
compared the simulation and HaloFit results for the con- 
vergence power spectrum (see Fig. 2 of the paper). The 
scale-dependence of the disagreement is qualitatively dif- 
ferent between real- and Fourier-space, reflecting the fact 
that the correlation function and the power spectrum is 
related to each other via the convolution. From vari- 
ous numerical tests we found that the shear correlation 
function remains accurate down to m ; n ~ 0.5 arcmin. 

Next we move on to the covariance matrix of shear cor- 
relation function. The covariance matrix describes how 
the correlation functions of different separation angles 
are correlate d with each oth e r. Fo llowi ng the methods 
develo ped in Uoachimi et al.l (|2008f ) and [Takada fc Jainl 
(2009), the covariance matrix for the correlation func- 
tion £+(6) is expressed as 

1 r°° 

Cov[£ + (0),£+(0')] = — - / ldlJo(W)J (l9')P K (l) 2 

1 poo poo 

+ THT m l'dl'J (W)J (l'6')T K (l,l'), (4) 

s J «/ 

where J7 S denotes the survey area and T K denotes the 
angle averaged lensing trispectrum. Note that we ig- 
nored the shot noise contribution due to random intrin- 
sic ellipticities. The derivation of Eq. (Q} assumes an 
ideal survey geometry; in other words Eq. @ is ap- 
proximately validated only for the case, 6,6' <C \/T^s 
and for large area surveys such as O s ^ 1000 deg 2 , as 



discussed in Appendix [X] in detail. Therefore it is ex- 
pected that for a small survey area the covariance ma- 
trix deviates from the formula above. We will call it 
the finite area effect whic h is examined in Appendix IA1 
See Uoachimi et all ()2008| ) for expressions of other covari- 
ance matrices Cov[£_(0),£_(0')] and Cov[£+(0),£_(0')], 
which are examined in § 13.21 

The first term in Eq. ((4]) describes the Gaussian contri- 
bution, while the second term is the non-Gaussian con- 
tribution. There is an important difference between the 
covariance of the convergence power spectrum and that 
of the real-space correlation function. Even for a pure 
Gaussian field, the first term in Eq. (U) is non-vanishing 
for the off-diagonal components with 6^6'. The corre- 
lation functions of different separations are always corre- 
lated with each other. Also note that the covariance does 
not depend on the bin width of angles. Thus, when us- 
ing the correlation function measurements for constrain- 
ing cosmological parameters, it is very important to have 
an accurate model of the covariance matrix in order to 
properly interpret the measurement. We will use Eq. (j4]) 
to compare the analytical prediction with the covariance 
measured from the simulations. 

2.2. Constructing a Gaussian Field from the 
Simulations 

The main goal of this paper is to quantify the relative 
importance of the non-Gaussian covariance to the Gaus- 
sian covariance as a function of angular scales and source 
redshifts. To address this, we first construct a Gaussian 
field using the ray-tracing simulations in order to sepa- 
rate out the Gaussian covariance contribution. We will 
hereafter call the constructed maps "the simulated Gaus- 
sian fields" . The reason why we use the simulated Gaus- 
sian fields instead of using the analytical prediction (the 
first term of Eq. |4} is as follows. Firstly, the ray-tracing 
simulations do not include large-scale modes beyond the 
simulation box size. Secondly, as we will show below in 
detail, the covariance measured from the Gaussian simu- 
lations shows a nontrivial dependence on the survey area 
that cannot be fully described by the first term of Eq. (|4j . 

We generated the Gaussian simulations according to 
the procedures below. Firstly, we Fourier-transformed 
each convergence field of 1000 ray-tracing simulations. 
Secondly, we make a Gaussian field by randomly select- 
ing each Fourier mode from the Fourier coefficients of 
1000 realizations, and then perform the inverse Fourier 
transform to obtain the real-space shear field imposing 
that the chosen Fourier modes satisfy the real number 
condition. Repeating this procedure, we made 1000 re- 
alizations of the Gaussian field. The simulated Gaussian 
fields generated in this way have the same power spec- 
trum on average as that of the original simulations, and 
contain Fourier modes over the same range of angular 
scales as in the original simulations. 

A justification of the Gaussian fields is given in Fig. [3J 
which shows the diagonal terms of power spectrum co- 
variance matrix, measured from the original and Gaus- 
sian simulated maps, as a function of multipoles. The 
simulation results are plotted relative to the expecta- 
tion for a Gaussian field, where the Gaussian covariance 
is equal to the squared power spectrum divided by the 
number of Fourier modes that are confined in each multi- 
pole bin. Therefore the deviations from unity arise from 



4 



Sato, Takada, Hamana and Matsubara 



- 1 1 1 1 1 1 1 1 1 1 1 II 1 

— 1 n 

z — 1 . u 

s 




* mode coupling 




■ no mode coupling 


X 


Gaussian 


- 


X 




X 








. « . ■ *• : a ■ ■ . ■ 




■ 

iM i 


1 ■ ■ 



10 s 10 3 
multipole 1 



Fig. 2. — The diagonal components of the convergence power 
spectrum covariance, divi ded by the expected Gaussian covariance 
(e.g. sec, Sato ct al. 2009), as a function of multipoles. The devia- 
tions from unity arise from the non-Gaussian errors. Note that the 
source redshift is z s = 1.0. The cross symbols are the simulation 
results from 1000 realizations, while the square symbols are the 
results obtained from the G auss ian simulations we generated from 
the 1000 simulations (see § 12.21 1. The Gaussian simulation results 
are consistent with unity over a range of multipoles. 

the non-Gaussian error contributions. The hgure explic- 
itly shows that, while the original simulations display 
stronger non-Gaussian covariances with increasing mul- 
tipoles, the Gaussian simulations are consistent with the 
Gaussian expectation over a range of multipoles we con- 
sider. 

Another justification is given in Fig. [3] The hgure 
shows the probability distributions of convergence k for 
one realization of the original ray-tracing simulations 
(red histogram) and the Gaussianized realization (blue 
histogram), respectively. The results are for source red- 
shift z s — 1.0. The original simulation, which includes 
the non-Gaussian contribution, shows a skewed distribu- 
tion, which strongly deviates from the Gaussian distri- 
bution with the same variance width. The distribution 
is better described by a log-normal distrib ution (dashed 
curve) , as studied in iTaruva et al.l (|2002f) , but shows a 
larger positive tail. On the other hand, the Gaussian 
simulation is in good agreement with the Gaussian dis- 
tribution (dotted curve). 

3. RESULTS: CALIBRATING THE NON-GAUSSIAN 
COVARIANCES 

3.1. Diagonal and Off-Diagonal Components of the 
Covariance Matrix 

In this section we study the covariance matrix of shear 
correlation function using 1000 ray-tracing simulations. 

The symbols in Fig.|4]show the diagonal components of 
the covariance matrix measured from 1000 simulations, 
for different source redshifts z s = 0.6,1.0 and 3.0, re- 
spectively. Again notice that each simulation realization 
has an area of 25 square degrees, so the plot shows the 
covariance expected when measuring the cosmic shear 
correlation from a survey with square-shaped survey ge- 
ometry and 25 square degrees. Here we ignored the in- 
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Fig. 3. — The probability distributions of convergence k com- 
puted from the original ray-tracing simulation and the Gaussian- 
ized simulation, respectively (taken from one realization of each 
simulations). The Gaussian simulation result is consistent with 
the Gaussian distribution that has zero mean and the same vari- 
ance as that measured from the simulation. 

trinsic ellipticity noise (we will come back to this later). 
First of all, we should stress that, by using the 1000 
realizations (25000 square degrees in total), we can ob- 
tain well-converged measurements of the covariance ele- 
ments over a range of the scales we cons ider. To be more 
quant itative, as shown in Appendix in iTakahashi et al.1 
(poll , covariances of the diagonal covariance elements 
scale approximately as (2/Ar) 1 / 2 , where -/V r is the num- 
ber of simulation realizations. Therefore the covariance 
elements are measured to 4%-level accuracies. 

We first compare the simulation results with the Gaus- 
sian error expectations (short-dashed curve) computed 
from the first term of Eq. ([4]). Note that we included 
only the Fourier modes confined within our simulations 
in the covariance calculation; we imposed P$(k) = at 
k < k{ as done in Fig. [1] Contrary to the result in Fig. [TJ 
the Gaussian predictions with the fc-cutoff effect appear 
to overestimate the simulation results on large separation 
angles, where the non-Gaussian errors are insignificant. 
On the other hands, the dotted curves show the results 
obtained from the Gaussian simulations we constructed 
(see § 12. 2|) , showing a good agreement with the simula- 
tion results on the large scales. Hence we conclude that 
the analytical prediction given by Eq. is not suffi- 
ciently accurate over the scales we consider. In fact, as 
studied in detail in Appendix[S] Eq. (|4]) can be valid only 
when the survey area is sufficiently large, such as 1000 
square degrees. In Appendix lAl we develop an empirical 
model to compute the Gaussian covariance taking into 
account the finite survey area effect. The long-dashed 
curve shows the modified analytical prediction for the 
Gaussian covariance, which is computed using Eq. (|A10[) 
as well as including the fc-cutoff for P$(k). The figure 
clearly shows that the modified analytical prediction well 
matches the Gaussian simulation results. Thus we need 
to properly account for both the finite range of Fourier 
modes and the finite survey area effect in order to obtain 
an accurate prediction of the covariance at large separa- 
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Fig. 4. — The diagonal components of the covariance matrix 
of shear correlation function for different source redshifts, z s = 
0.6, 1.0 and 3.0. The cross, square and circle symbols are the simu- 
lation results from 1000 realizations for each source redshifts. The 
solid curves show the halo model predict ions including th e 1-halo 
term and the halo sample variance (see ISato et~aT] [20091, for the 
details). The dotted cu rves arc the results obtained from the Gaus- 
sian simulations (see § 12.21 and Fig. [2j. The short-dashed curves 
show the analytical predictions computed from the first term of 
Eq. i4l . where we included only the Fourier modes contained within 
our simulation boxes, i.e. imposed the condition P$(k) = at 
k < kf in the covariance calculation as in Fig. [l] On the other 
hand, the long-dashed curves are the results which are computed 
using an empirical model to account for the finite survey area effect 
(see Appendix [X] for details) . This modified Gaussian prediction 
better matches the Gaussian simulation results. 



tions. 

The simulation results (symbols) start to deviate from 
the Gaussian simulation results on small separation an- 
gles due to the non-Gaussian error contribution. The 
plot shows that the non-Gaussian errors are more sig- 
nificant on smaller angles and for lower source redshifts, 
because of stronger nonlinear clustering in structure for- 
mation at lower redshifts and on smaller length scales. 
For comparison, the solid curves show the halo model 
predi ction developed in our previous paper (|Sato et al.l 
2009). To be more precise, the halo model is used to 
compute the non-Gaussian covariance for the assumed 
cosmological model, inc l uding the "halo sampling vari- 
ance term" (jSato et al.l 120091 ) in addition to the non- 
Gaussian term (1-halo term). Then the solid curves are 
the sum of the halo-model-computed non-Gaussian co- 
variance and the Gaussian covariance that is computed 
from the Gaussian simulations. The halo model predic- 
tions qualitatively well reproduce the simulation results 
on small separations as well as the source-redshift de- 
pendence. However, the halo model also shows a sizable 
disagreement due to the limitation. 

Fig. [5] shows the covariance matrices comparing the re- 
sults for the original ray-tracing simulations, the Gaus- 
sian simulations and the Gaussian prediction computed 
from the first term of Eq. (fj|, where the fc-cutoff and 
the finite survey area effect are included. Note that all 
the results are for z s = 1. The Gaussian prediction 
(right panel) and the Gaussian simulation results show 
a nice agreement for the scale-dependences and ampli- 



tudes. The original simulation results show totally dif- 
ferent scale-dependences from the Gaussian results and 
display greater amplitudes at smaller separation angles 
(left-lower corner) due to stronger non-Gaussian contri- 
butions. 

Next we study the correlation coefficients to quantify 
strengths of the off-diagonal covariance components rel- 
ative to the diagonal components: 



Cov[£ + (0),e+(0')] 



vcovK+(0),M0)]cov[Mne4-(0')] 

(5) 

The correlation coefficient is defined so that r — 1 for 
the diagonal components when 9 — 9' . For off-diagonal 
components when 9 ^ 9' , r — > 1 implies strong corre- 
lation between the two correlation functions of different 
angular scales, while r — >■ means no correlation. 

Fig. [6] shows the correlation coefficients 
r[£+(0), £+(#')]. There are significant correlations 
between different separation angles. Comparing the 
dotted and solid curves manifests less significant cross- 
correlations in the non- Gaussian errors than in the 
Gaussian errors, implying that the non-Gaussian errors 
preferentially contribute to the diagonal components. 

3.2. A Fitting Formula of the Non-Gaussian 
Covariance Elements 

By using the simulation results shown up to the pre- 
ceding section we derive a fitting formula to compute 
the non-Gaussian covariance as a function of separation 
angles and source red shifts. To do th i s we e mploy the 
method developed in iSemboloni et al.l (|2007) , which is 
to derive the calibration function that gives the non- 
Gaussian covariance contribution relative to the Gaus- 
sian covariance. O ur results bring several i mprovements 
over the results in Se mboloni et al.l (|2007l ). Firstly, we 
use 1000 ray-tracing realizations, for each source redshift 
[zs = 0.6, 0.8, 1.0, 1.5, 2.0, 3.0). Secondly, we carefully es- 
timated the Gaussian covariance contribution using the 
Gaussian simulations (see § I2.2|) . thereby enabling us to 
reliably estimate the relative contribution of the non- 
Gaussian covariance. 

We define F(9,9';z s ), the ratio of the non-Gaussian 
covariance relative to the Gaussian expectation for the 
covariance matrix: 



F(8,6';z s ) 



Cov NG K + (flU + (g');^ 
Cov G K+(0),£f(0');*.] 



(6) 



Here we used the 1000 simulations to compute the Gaus- 
sian and non-Gaussian covariance matrices appearing in 
the numerator and denominator of the above equation. 
As shown in Figs.|4]and[5]the non-Gaussian contributions 
are important only on small separation angles, 10 ar- 
cmin, for a ACDM cosmology and for source redshifts we 
consider. 

Fig. [7] shows the simulation results for the diagonal 
components of F(9, 9') for different source redshifts. The 
fitting formula for F(9,9') is given in Appendix iBl The 
figure shows that all the curves go to unity, F(9,9) = 1, 
on very large separations as expected; the non-Gaussian 
covariances become negligible on such l arge separations. 
This c an be contrasted to the results in ISemboloni et all 
(|200l (see Fig. 1 in their paper), where the correspond- 
ing curves do not go to unity, even go below unity at 
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1 10 100 1 10 100 1 10 100 

Ofarcmin) 0(arcmin) O(arcmin) 



Fig. 5. — The covariance matrices Cov[£-|_ (8), (6')] as a function of separation angles 8 and 8' , for the ray-tracing simulation realization 
(left panel), the Gaussian simulation (middle) and the analytical prediction computed from the first term of Eq. JJJ (right), respectively. 
All the results are for source redshift z s = 1.0. For the analytical Gaussian covariance, we included the finite survey area effect as well as 
the fc-cutoff in the power spectrum calculation, as in Fig. [4] 




Fig. 6. — The correlation coefficients r[£_|-(0), £-)_((?')] as a function of separation angles 8, where 8' is kept fixed to 8' = 0.794, 7.94 and 
79.4 (arcmin) in the left, middle and right panels, respectively. In each panel we show the results for different redshifts, z s = 0.6, 1.0 and 
3.0. The solid curves denote the halo model predictions, while the dotted curves denote the Gaussian error predictions. 



large separations. This is because we carefully com- 
puted tlve_GaussjaiiTOvariance contribution (see § I2.2[) . 
while fSembolo ni et al.l (|2007T ) used the analytical predic- 
tion (the first term of Eq. [j) to estimate the Gaussian 
covariance, which turns out t o overestimate the sim ula- 
tion results. Also note that in lSemboloni et all (|2007l ) all 
the Fourier modes setting the lower bound to k = or 
/ = are included in the Gaussian covariance calculation 
(Eq.SD. 

The figure also shows that the non-Gaussian errors are 
significant on smaller separations and for lower source 
redshifts, where nonlinear clustering is more evolving. 
For a source redshift z s ~ 1, which is a typical depth of 
the Subaru- type survey, F £ 10 on scales smaller than 
1 arcminutes, meaning that the non-Gaussian contribu- 
tion to the diagonal covariance is greater than the Gaus- 
sian contribution by more than a factor of 10. The 
non-Gaussian covariance amplitudes are by accide nt sim- 
ilar to what is found in iSemboloni et al.1 (j2007j ). even 
though the power spectrum normalization as is quite 
di fferent: as = 0-7 6 in o ur simulations, while as = 1 
in ISemboloni et al.l (120071). Ag a in thi s is subscribed to 
the fact that ISemboloni et al.l (|2007l ) under-estimated 



the non-Gaussian error contribution. The normalization 
as = 0.76 we assumed is slightly lower than the cur- 
rently most-favored value, a s = 0.8 (see WMAP 7- year 
result. iKomatsu et al.|[201lT ) . Therefore the fitting func- 
tion F(8,9') we calibrated may slightly underestimate 
the non-Gaussian errors (by about 20%). To correct for 
this, we can use the halo model prediction to estimate 
the difference in the non-Gaussian covariance amplitudes 
for different as values, and then multiply the correction 
factor with F(0,0'). 

Let us summarize how to obtain the covariance matrix 
for given survey parameters from our fitting formula: 

1. Compute the first term of Eq. (|Al0j) to estimate 
the Gaussian covariance matrix, for given survey 
area and source redshift. 

2. Use the fitting function Eq. JET) for F(0,9'; z s ) to 
obtain the non-Gaussian covariance contribution. 

3. Multiply the quantities in the steps 1 and 2 to ob- 
tain the total covariance matrix. 

Exactly speaking, the survey-area dependence of the non- 
Gaussian covariance is not as naively expected; Cov oc 
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Fig. 7. — The diagonal elements of shear correlation covariancc 
relative to the Gaussian covariance contribution, as a function of 
separation angles and source redshifts (see Eq. [6] for the defini- 
tion). The results arc measured from the 1000 realizations, where 
the Gaussian contributions arc estimated from the 1000 Gaussian 
simulations as in Fig. [2] 



l/fl s does no t hold due to th e halo sampling variance 
contribution ()Sato et al.l I2009T > . However we have found 
that the residual dependence is relatively small, and our 
fitting formula for the covariance matrix is approximately 
valid for a survey area we are most interested in ( ^ 100 
deg 2 ) ; 

While we have so far ignored the intrinsic ellipticity 
contribution to the covariance, we need to further in- 
clude the effect. First of all, the intrinsic ellipticity noise 
contributes only to the Gaussian covariance, and does 
not affect the non-Gaussian covariance, as long as the 
intrinsic ellipticity alignment can be ignored. Hence we 
can easily include the shot noise contribution as follows. 

One method is a simulation based method. We can 
generate the Gaussian field including the shot noise con- 
tribution by replacing the power spectrum P K (l) with 
Pk(1) + of/^g m generating the simulation field, where 
<7 e is the rms intrinsic ellipticity per component and n g 
is the mean number density of source galaxies. Here 
we need to assume a e and n g that we measure from 
galaxies in a given survey. If we measure the covari- 
ance matrix from the Gaussian simulations generated in 
this way, the covariance matrix includes the contribu- 
tion arising from the term oc (a 2 /h g ) 2 as well as the 
mixed term between the cosmological Gaussian field and 
the shot noise term arising from the term o c P K (l)a 2 /n g 
(jSchneider et al.ll2002t Ijoachimi et al.ll2008ri. 

Ano ther way i s following the method in ljoachimi et ail 
(120081) (also see iSchneider et al.|[2002h . which developed 
an analytical formula of the Gaussian covariances includ- 
ing the intrinsic ellipticity noise contribution. To be more 
explicit, the formula is given in terms of the convergence 
power spectrum P K (l), as in the first term of Eq. (j4]). 
Then the intrinsic noise contribution can be incorporated 
by replacing P K (l) 2 with (P K (l) +a 2 /n g ) 2 . Here we need 
to assume er e and n g that are measured from a given 
survey. 
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Fig. 8. — The diagonal components of shear correlation covari- 
ances Cov[£_(0),£_(0')] and Cov[£+(0), as a function of 
separation angles at source redshift z s = 1.0. The cross symbols 
are the simulation results measured from the 1000 realizations. The 
dotted cur ves a re the results obtained from the Gaussian simula- 
tions (see § 12.21) . while the dashed cur ves the theoretical pre dictions 
computed from Eq. (30) or Eq. (31) in Joachimi et al. (2008), where 
we included the fc-cutoff in the lensing power spectrum calculation. 



Finally we comment on other contributions to the 
covariance matrix we have so far ignored, which 
are the covariance contributions arising from another 
shear correlation function £-(#): Cov[£_(0),£_(0')] and 
Cov[£+(0),£-(0')]. Fig. shows the results. Again 
the non-Gaussian error contributions are significant 
at scales ,$ 10'. Compared to Fig. 21 the covari- 
ancc matrices Covj£_(0),£_(0')] and Cov[£ + (0),£_(0')] 
have smaller amplitudes than Cov[£+(0),£+(0')] does: 
for example, the amplitudes of Cov[£_(0), £_(0')] 
and Cov[£ + (0),£_(6>')] are smaller than that of 
Cov[£ + (0), £+(#')] by a factor of 100 and 10, respectively. 
However, the genuine effects need to be understood in 
terms of the signal-to-noise ratios that are roughly es- 
timated as £+/Cov[£ + , £+] or £ + £_/Cov[£+, £_] at each 
separation angles. Since the correlation function £_ has 
smaller amplitudes than £+ does, by a factor of 10 at 
separ ations ~ 1' (e.g. see Fig. 2 in Schneide r et al.1 
2002), therefore the covariances Cov[£_(0),£_(0')] and 
Cov[£+(0),£_(0')] are not negligible. 

We tried to derive fitting formulas for the non- 
Gaussian covariance contributions to Cov[£_(0),£_(0')] 
and Cov[£+(0),£_(0')] in the similar manner as done in 
Eq. ((6]). However, due to complex scale-dependences of 
the Gaussian covariances as implied in Fig. [SJ we could 
not find useful fitting formulas that are expressed by sim- 
ple analytical functions. Therefore, the covariance ma- 
trix contributions need to be directly calibrated from the 
ray-tracing simulations. The table-format covariances 
for Cov[£_(0),£_(0')] and Cov[C+(0),f_(0')] are avail- 
able upon request. 

4. CONCLUSION 

In this paper, we have developed the theoretical model 
of the covariance matrix of cosmic shear two-point cor- 
relation function taking into account the effect of finite 
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survey area and the effects of non-linear gravitational 
clustering in large-scale structure. 

We found that the survey-area dependence of the 
Gaussian covariance, which scales as oc l/fi s in the com- 
monly used formula, does not hold for a small area sur- 
vey due to the effect of finite survey area. The conven- 
tional formula is valid only when the survey area is suffi- 
ciently wide such as 1000 square degrees. We examined 
the residual s urvey-area depe n dence using the method 
developed by iSchneider et al.l ((2002), and obtained an 
empirical formula which reproduces our results. 

We examined the non-Gaussian covariance as a func- 
tion of angular scales and source redshifts by using two 
sets of simulation data: (1) the ray-tracing simulations 
for the standard ACDM cosmology and (2) Gaussian 
fields which have the same power spectrum on average as 
that of the original simulations. The non-Gaussian errors 
become more significant on smaller scales and at lower 
redshifts. We compared the simulation results with halo 
model predictions and found that the halo model qual- 
itatively well reproduces the non-Gaussian error over a 
wide range of separation angles and for redshifts we have 
considered. However, the halo model also displays sizable 
disagreeme nt with the simulation r esults. 

Following iSemboloni et al.l (|2007l) , we derived the cal- 
ibration function to compute the non-Gaussian covari- 
ance contribution relative to the Gaussian covariance 
(Eq. [6|). The Gaussian field data allows us to cleanly 
separate out the Gaussian covariance contribution accu- 
rately, thereby enabling us to reliably estimate the cal- 
ibration function. We found that the calibration factor 
at arcminute scales can be high as ~ 20 and ~ 10 for 
source redshifts of z s = 0.6 and 1.0, respectively. The 



transition between Gaussian and non-Gaussian covari- 
ance occurs around 10 and 5 arcminute for z s —0.6 and 
1.0, respectively. Therefore, when one derives constrains 
on the cosmological parameters from cosmic shear cor- 
relation measurements, it is important to properly ac- 
coun t for the non-G auss ian effects. The fitting formulae, 
Eq. (IA10p and Eq. (IB1[) . developed in this paper allow 
one to compute the covariance matrix including the non- 
Gaussian contribution for given survey parameters (the 
survey area and the source redshift). 

Simulation data (1000 convergence power spec- 
tra and cosmic shear correlation functions for £+(9) 
and £_(#).) are available u pon request (contact 
masanori@a. phys.nagoya-u.ac.jp). 
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APPENDIX 

A. FINITE AREA EFFECT OF THE GAUSSIAN COVARIANCE 

In this section, we study the validity of the first term of Eq. ((4]), which gives the Gaussian error prediction for the 
shear correlation function covariance, and will show the prediction is only valid for a large-area survey covering more 
than 1000 square degrees. For this purpose we will use the method developed in ISchneider et al] (|2002l ). 

Let us begin with considering an estimator of the shear correlation function of a separation angle 9, £+(#)• For a 
given galaxy catalog the shear correlation function can be estimated as 

£+(0) = -T?— EE ( e *(«) +e x{l) e xU) )A(\9 t -6 3 \;9), (Al) 

J* pair 

% j 

where we have used the abbreviated notations such as e t u\ to denote the ellipticity component of the i-th galaxy at 
the angular position 8i, iVpair is the total number of pairs of galaxies that are separated by the separation angle 9 
and the index i in summation runs over all the galaxies in the catalog. The tangential and cross components of the 
ellipticity e = t\ + \t2 at position # are defined as 

e t = -Ke(ee~ 2 ^); e x = -1m (e e" 2 ^) , (A2) 

where IZe and Xm denote the real- and imaginary-parts of the quantities and ip is the polar angle of the separation vector 
between two galaxies, 6. The component, et, is defined as the ellipticity component in parallel or perpendicular direction 
relative to the line connecting the two galaxies. On the other hand e x are measured from the ±45° rotated components. 
The function A(|0< - 8j\;d) is a selection function defined in that A(|0, - 8j\;9) = 1 if 9- 59/2 < \8 t -8j\< 9 + 59/2, 
otherwise A(|0j — 8j\; 9) = 0, where 59 is the bin width of separation angle. The total number of galaxy pairs is given 
as A^pair = £\. A(\8i — 8j\;9). The ensemble average of the estimator Eq. (|A1|) is found to indeed give the shear 
correlation function: 

(i+m = ^Y,( £ ^tU)+^ ) e xU) )A(\d l -8 ] \-,9) = J—^t + {\e i -e J \)AQO i -O j \-,0)=t + (6). (A3) 

^*pair . . ^'pair ■ . 
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Similarly the covariance is denned in terms of the estimator £ + (9) as 

Cov[Z + (0te + (O f )] = (£ + (e)i + (9')) Z+m+W). (A4) 

For simplicity, let us co nsider the diagonal pa rts of the covariance matrix, 6 = 6'. For the Gaussian field, the diagonal 
parts are computed as (|Schneider et al.ll2002t ): 

Cov[£+(0U + (0)] = ^-E<( e '«%-) +e =<« £i <0-)) (e t(i) e t(m) + e><(oe><(m) ))A(|0 4 -^|;0)A(|^-0 m |;0)-e+(0) 2 

pair ij[ m 

= W7 EK+d^ - d m\)Z+{\0 3 - 0i\) + cos[4(^ rn - tpuM-QOi - 6 m \)t-(\0i - Oi\W\Oi - e j \ ] 9)A(\e l - e n \-,6) 

pair ijl m 

= J^2 E Fijlmi (A5) 
P air ijlm 

where ifi m is the polar angle of the difference vector 8i — 8 m , and we used the fact that an estimator for £-(0) is 
analogously defined as 

L(0) = ^— E ( e *W e *0) ~ e x W ex )) A(|0 t -0 3 \-9). (A6) 

^*pair 

Now we use Eq. (|A5|) to study the effect of finite survey area. This can be done by comparing the result with the 
first term of Eq. Q because Eq. ^ ignores the survey geometry effect. To do this we performed the simplified test. 
Firstly, we randomly distribute galaxy positions with square shape geometry. Then, rather than working on the shear 
field in simulations, we will compute the summation in Eq. (|A5I) by using the tabulated data of £+(0) and £-(#). We 
estimate £+(#) and using Eq. (TT]) and Eq. ([2]). In this case, the summation such as runs over all the galaxies 

in the square-shaped mock simulation, and the separation angle such as \0i — 0j\ can be exactly computed from the 
separation between the two galaxies chosen. 

The symbols in Fig. IA1I show the results of Eq. (|A5[) for the diagonal covariance matrix elements as a function of 
separation angles and mock simulation areas. The solid curves arc the results computed from the first term of Eq. Q. 
In both cases we used the HaloFit prediction to compute the lensing power spectrum, which is then used to compute 
the shear correlation function as well as to compute the first term of Eq. (Q} . Note that we do not plot the symbols 
where separation angle becomes comparable with the scale of the mock simulation area due to avoid the boundary 
effect of the square-shaped mock simulation. 

The figure clearly shows that the symbol is in good agreement with the solid curve for the largest area we consider, 
f2 s = 3200 deg 2 , however, the two results disagree in amplitudes for the smaller areas, although the shape looks 
similar. In addition, comparing the results for O s = 1.56 and 25 deg 2 manifests that the disagreement is greater for the 
smaller survey area. This disagreement arises due to the effect of the finite survey area. Hence the conventionally used 
expression for the Gaussian covariance, the first term of Eq. (jj), appears to overestimate the covariance amplitude. 
The formula is only valid for a sufficiently large area such as il s 1000 deg 2 . In other words, due to the finite survey 
area effect, the covariance of shear correlation function does not scale with survey area as Cov oc l/f2 s , which is 
assumed in the conventional formula (Eq. [4]) . 

From now on, we will in more detail study how the difference between the results of Eqs. (Q| and (|A5|) ari ses. Let 
us consider two survey areas A and B, assuming the area B is included in the area A (A 3 B). By using Eq. (|A5|) the 
diagonal covariance matrix for a survey A can be expressed as 



Cov A [£ + (0U + (0)] = a . 



air \ijlmeA 



^ ^ Fijlm ~l~ 4 ^ ^ ^ijirn ~l~ 2 ^ ^ -^V/'/m 

ijimeA-B ieB ij&B 

j7meA-B ; mS A-B 



i/GB imEEB j/m(EB ijlm^'B 

jmeA-B jieA-B teA-B / 



(A7) 



where the notation ij 6 B is introduced to mean that the indices i and j are included in B, and the notation Im € A — B 
means that the indices I and m are included in the region that is in A, but not in B. We will consider the case that 
survey area A is sufficiently large in the following discussion. We define the value arising from finite area effect as 

FAE = gCovA [£+(#) ,£+(#)] — Covb [£+(#)> £+(#)]• (A8) 

If B is close to A, there is no fini te a rea effect and FAE should be zero because of survey area as Cov oc 1 / O s . If B is 
close to A, the final term of Eq. (|A7[) has a dominant contribution because the number of galaxy pairs is greater than 
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Fig. Al. — A numerical test for studying the finite area effect on the Gaussian covariance of sh ear correlation (see text for the details). 
The symbols show the results for the diagonal covariance elements, which are computed using Eq. I IA5I I. for different areas of square-shaped 
mock simulations: U s = 1.5625 deg 2 , 25 deg 2 and 3200 deg 2 . The source redshift is z s = 1.0. The solid curves are the prediction computed 
from the first term of Eq. (|4} . The two results agree for a large area such as Q B = 3200 deg 2 , but disagree for the smaller areas probably 
due to the finite-area effect. 



other terms. In this case Eq. (|A8|) becomes 



1 



FAE = £cov A [£ f (d),e f (0)] - COVBE+ (*),£+ (*)] * ^TF^ £ F ^ ~ Covb[£+(0),£ + (0)] 



i 



-^pair ly/mel 



F ij i m -Cov B [Z + {0),{+(d)] = 0. 



(A9) 



Thus Eq. (|A8j) has an asymptotic limit of Cova = Covb when A = B. However, the other terms in Eq. (|A7j) are not 
negligible if A ^ B, and the covariances for a general case do not scale with survey area as Cov oc l/f2 s , due to a finite 
survey area effect. 

Therefore we will estimate a fitting formula to account for the survey area dependence by assuming the form 



Cov[£+(0), £+((?')] 



1 



1 



Idl / l'dl'J (W)Jo(l'6')T K (l,l'), 

I! JO 

(A10) 



where f(£l s ) denotes the new survey area dependence. We parametrize /(f2 s ) as 

oc(z a ) 



/(Os) 



(All) 



We estimate the 2 parameters using the mock simulations as in Fig. IA11 for different source redshifts z s 
0.6, 0.8, 1.0, 1.5, 2.0 and 3.0. The mock simulation results are well fitted by the parameters 



a(z s ) = a ± z c s 



(3(z s ) = fazf* 



(A12) 



The best-fitting parameters are found to be (ai,a 2 ) = (3.2952,-0.316369) and (/?i,^ 2 ) = (0.170708,-0.349913), 
respectively. As shown in Fig. IA2l we made the fitting over survey areas of 1.5625 < f2 s < 1600 deg 2 using the diagonal 
covariance at 50 arcmin. However we checked that the above fitting formula fairly well reproduce the results for 
different separations and for the off-diagonal covariance elements. A caution on the use of the fitting formula is the 
output value should be replaced to unity if the value is below unity. 

B. FITTING FORMULA FOR CALIBRATION MATRIX 

In this section we estimate the calibration function F(Q, 9; z s ) w e stud ied in § 13.21 which can be used to estimate the 
non-Gaussian covariance in combination with the fitting formula (|A10I) for the Gaussian covariance. We parametrize 
F(9, 9'; z s ) as 

F(9,9';z s )^(a(z s ) + -^-) ■ d(z s f-°X (Bl) 



?')C(2.) 

The 4 parameters are estimated from the ray-t raci ng simulations for 6 different source redshifts z s = 0.6, 0.8, 1.0, 1.5, 2.0 
and 3.0. As demonstrated in Figs. IB1I and IB21 the simulation results are well fitted by the following best-fitting 
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Fig. A2. — A fitting formula for the survey area depen denc e of the Gaussian covariance, /(O s ) in Eq. tAlOt . The sy mbols are the mock 
simulation results for different source redshifts as in Fig. lAll and the solid curves are the fitting formula given by Eq. jAlU . We did the 
fitting by matching to the mock simulation results for the diagonal Gaussian matrix at the scale of 50 arcmin. 




0(arcmin) 

Fig. Bl. — The diagonal components of the calibration function to compute the non-Gaussian covariance, F(8 , 
definition). The symbols are the same as in Fig. [3 and the solid curves are the fitting formula given by Eq. |)B1J| , 



,) (see Eq. [6] for the 



parameters: 



a(z s ) = -z" 1 exp(a 2 z s ) + a 3 
b(z s ) = b lZ b / + h 
c{z s ) = cizf + c 3 
d{z s ) = d lZ ^ + d 3 , 



(B2) 



where (01,03,03) = (-3.7683,0.9752,1.4048), (&i,& 2 ,& 3 ) = (10.7926,-2.0284,-0.2266), (ci,c 2 ,c 3 ) = 
(-0.3664,-0.5733,0.6863), and (d 1: d 2 ,d 3 ) = (0.2450,0.1218,0.7076), respectively. We made this fitting over an- 
gular scales 0.5 < {9, 0'} < 10 arcmin. Figs. IB 1 1 and IB2I clearly show that the fitting formula above well reproduce the 
simulation results, and the accuracy of the fitting formula is within about 25%. However it should be noted that the 
fitting formula is only applied to source redshift ranges 0.6 < z s < 3.0 and angular scales below 10 arcmin. This is 
sufficient because the non-Gaussian covariance is important only below 10 arcmin. 
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